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We introduce a real time version of the functional renormalization group which allows to study 
correlation effects on nonequilibrium transport through quantum dots. Our method is equally 
capable to address (i) the relaxation out of a nonequilibrium initial state into a (potentially) steady 
state driven by a bias voltage and (ii) the dynamics governed by an explicitly time-dependent 
Hamiltonian. All time regimes from transient to asymptotic can be tackled; the only approximation 
is the consistent truncation of the flow equations at a given order. As an application we investigate 
the relaxation dynamics of the interacting resonant level model which describes a fermionic quantum 
dot dominated by charge fluctuations. Moreover, we study decoherence and relaxation phenomena 
within the ohmic spin-boson model by mapping the latter to the interacting resonant level model. 
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I. INTRODUCTION 

Obtaining a deeper understanding of nonequilibrium 
phenomena in the presence of many-body correlations is a 
major challenge in condensed matter physics. One partic- 
ularly well-defined and controllable setup - both from an 
experimental and a theoretical perspective - are so-called 
quantum dots which feature a few correlated electronic 
degrees of freedom coupled to noninteracting leads. Their 
low-energy equilibrium physics is typically governed by 
the appearance of an energy scale which compared to the 
bare scales is strongly renormalized through the Coulomb 
interaction. A prominent example is the Kondo efFect>i 
If a (nearly) odd number of electrons reside on the dot, 
spin fluctuations are strong below the Kondo tempera- 
ture Tk- The latter depends exponentially on the bare 
system parameters given by the local charging energy 
U, the level- lead hybridization F ('kinetic energy'), and 
the level position e. Another example are quantum dots 
dominated by correlated charge fluctuations, in which the 
decay rate F is renormalized. A prototypical model to 
describe such a scenario is the interacting resonant level 
model (IRLM),^ 

In equilibrium, methods based on the renormalization 
group (RG) idea proved to be powerful tools to treat 
many-body quantum dot problems with a hierarchy of en- 
ergy scales;^"— Several RG(-like) approaches were there- 
fore recently extended in order to address three different 
nonequilibrium scenarios. In the most simple case one is 
interested in the nonequilibrium steady state induced by 
coupling to two (or more) leads that are held at chemical 
potentials differing by a bias voltage V. The Hamilto- 
nian itself is taken as time-independent. Prominent ob- 
servables are the current J through and the occupancy 
fi of the dot levels. Problems of this class were stud- 
ied using analytical RG(-like) methods such as the poor 
man's RG,^^ the real-time RG (RTRG),^^ and the flow- 
equation approach.— Numerical methods include Wilson's 



numerical RG (NRG) frameworkiSiii (note that it was 
recently questionec^^ if NRG's inherent logarithmic dis- 
cretization is reasonable in nonequilibrium), the time- 
dependent density-matrix renormalization group,— an 
iterated path integral approach,— and time-dependent 
quantum Monte Carlo^ (the last three are not RG- 
based). 

A more intriguing (but certainly more complex) task is 
to study the relaxation dynamics towards a steady state 
configuration, i.e. to ask: How does a specific nonequi- 
librium state time-evolve under a still time-independent 
Hamiltonian? RG(-like) methods employed to ad- 
dress this question are the flow-equation approach,— 
RTRG,.^'i^ and NRG.iSiii The last - and for a theoreti- 
cal description yet more challenging - class of problems 
are those in which the Hamiltonian carries an explicit 
time dependence; in the context of quantum dots, charge 
pumping is a typical examplej^^ The RTRG was recently 
extended to investigate this scenario.— 

The functional renormalization group (FRG)^ imple- 
ments Wilson's RG idea in terms of an a priori exact 
infinite hierarchy of differential flow equations for the 
many-body vertex functions. It has distinct general ad- 
vantages over other Wilson-like RG procedures: Func- 
tional RG (i) can be applied directly to microscopic mod- 
els and not only to effective field theories, (ii) provides 
information on all energy scales and not solely on the 
low-energy limit, and (iii) allows for a flexible introduc- 
tion of the flow parameter (cutoff). The key approxima- 
tion is to truncate the infinite hierarchy at a given order. 
This is controlled for weak to intermediate interactions 
- a parameter regime which for certain problems might 
still be dominated by electronic correlations. Indeed, the 
FRG was shown to provide a reliable tool to study the 
linear-response physics of single- and multi-level quan- 
tum dot setups.—"— It was recently extende d^^'^^ in or- 
der to investigate the steady-state limit in nonequilib- 
rium (i.e., to treat the first class of problems discussed 
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above) where it captures certain aspects of nonequilib- 
riuni Kondo physics^^ and yields a comprehensive pic- 
ture of the finite-bias transport through a quantum dot 
dominated by correlated charge fluctuations.— The sin- 
gle impurity Anderson model as well as the IRLM were 
employed as prototypical examples. 

This paper aims at a natural but nontrivial gen- 
eralization of the functional RG framework which al- 
lows to tackle real time relaxation dynamics as well as 
time- dependent Hamiltonians of interacting quantum dot 
problems. In complete analogy with the prior extension 
of the FRG from linear response to steady-state nonequi- 
librium, this requires (i) to derive the exact hierarchy of 
flow equations from a functional that generates real-time 
vertex functions, (ii) to introduce a cutoff which preserves 
symmetries (such as causality) in nonequilibrium, (iii) to 
consistently truncate the infinite hierarchy at a given or- 
der and to formulate a closed set of flow equation, and 
finally (iv) to implement an algorithm which solves them 
'numerically exact' in reasonable time. In contrast to all 
prior applications, devising such an algorithm is involved 
for the problem at hand. To this end, we organize our 
paper as follows: 

In Seciniwe introduce a general Hamiltonian that de- 
scribes an at this point unspecified quantum dot tunnel- 
coupled to noninteracting leads. We allow for time- 
dependent dot and tunnel parameters. In order to even- 
tually set up the FRG flow equations, we discuss some 
basics of nonequilibrium single-particle Keldysh Green 
functions that depend on two time arguments. We 
shortly illustrate how to express the current J and dot 
occupancy n in terms of those quantities. 

In a next step (Sec. IIIip we derive the exact hierarchy 
of flow equations and discuss their form after truncating 
at the lowest nontrivial level. Even though this approx- 
imation - which is the only approximation within our 
approach - can be strictly motivated for small Coulomb 
interactions only, it was successfully used to describe 
aspects of correlation physics for a variety of quantum 
dot setups in equilibriuroSi as well as in steady-state 
nonequilibriumi^i2i We stress that the details of our 
FRG implementation for fermions are different from ear- 
lier extensiona^ii^ of the method to study the time evolu- 
tion of interacting bosons. Importantly, we do not rely on 
the generalized Kadanoff-Baym ansatz to approximately 
solve the Dyson equation (which was used in Ref. l2a ) 
since wc are sceptical about its general validity. 

In Sec. IIVI we specify our (so far general) FRG ap- 
proach for the interacting resonant level model. The lat- 
ter describes a single spinless fcrmionic level with energy 
e{t) that is locally coupled to two Fermi liquid leads via 
a tunnel matrix element t(<) and a Coulomb interac- 
tion U{t). In the steady-state limit, the corresponding 
flow equations could partially be solved analytically (in 
frequency space) .— This is no longer possible in our real- 
time representation, and we need to resort to a numerical 
treatment. Implementing a 'numerically exact' solution 
efficiently on a standard computer is not straightforward. 



We present aspects of our algorithm in Sec. |Vl 

Section I VII is devoted to elaborating the current J{t) 
and the dot occupancy n{t) of the IRLM obtained from 
our FRG scheme. We focus on time-independent (but 
otherwise general) system parameters for reasons of sim- 
plicity but emphasize that the framework directly allows 
to study an explicit time dependence (results will be pub- 
lished elsewhere^). In a nutshell, wc find two different 
rcnormalized relaxation rates, characteristic oscillations 
with frequencies given by the level position relative to 
the left and right lead chemical potentials, and power-law 
corrections to the exponential time dependences. We do 
not observe the appearance of secular terms frequently 
encountered within perturbation theory.— A comparison 
of our results to real-time RG dat a^^i'^^ strongly supports 
that the lowest-order FRG fiow equations in real time 
capture hallmarks of correlation physics within the in- 
teracting resonant level model. 

In absence of a bias voltage, the IRLM can be mapped 
onto the ohmic spin-boson model^ whose relaxation 
dynamics we investigate as a second application (Sec. 
IVII|) . We relate our results to predictions from field 
theory^! and the so-called improved noninteracting blip 
approximation.— 

A summary and a perspective for future applications 
of our time-dependent functional RG scheme are given 
in Sec. IVIIII In the Appendix we outline some technical 
details of the calculation of the Keldysh Green functions. 

II. OPEN FERMI SYSTEMS IN THE KELDYSH 
FORMALISM 

Hamiltonian — We aim at discussing time-dependent 
nonequilibrium transport through a quantum dot cou- 
pled to two or more leads (i.e., transport through an open 
Fermi system). To this end, we introduce the Hamilto- 
nian 

Hit) = i/^°'(o + ^[i/r + Hrnm ■ (i) 

a 

The dot part constitutes of a single-particle term and a 
two-particle interaction, 

H'^°\t) = iJo^°'(i) + H"'\t), (2) 
H^''\t) = Y,e,,it)dld„ (3) 

H'^\t) = Wu,,ki{t)4dl'^idk, (4) 
4 ijkl 

where we employ standard second quantized notation. 
The reservoirs a are modelled as noninteracting, 

Hr = T^kA^ck^^ (5) 

and they arc tunnel-coupled to the dot through 

H^^t)=Y,JkAt)cld^ + ll.c.. (6) 

ka ,2 
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We explicitly allow for a time dependence of the param- 
eters of i/'^"* and i?™"P. 

Initial statistics — We assume that the system is pre- 
pared using a product density matrix p at time t = 
- a situation which arises naturally when the dot and 
the reservoirs are decoupled for t < 0. Furthermore, the 
reservoirs are supposed to initially be in grand canonical 
equilibrium with temperature Tq, and chemical potential 



p(t = 0) = po = pr ® pZo ® • ■ • ® pZ,o, (7) 



res 
Pa.O ■ 



where Nn 



We choose units with fee = 1, 



h = \, and electron charge e = 1. Finally, we assume 
that the statistical operator Pq°* at i = allows for the 
application of Wick's theorem^ and that its matrix rep- 
resentation commutes with the initial single-particle dot 
Hamiltonian: [pg°*, e(t = 0)] = 0, where e is the matrix 
with entries . For an initially empty quantum dot these 
requirements are trivially fulfilled; it is then irrelevant 
whether or not the two-particle interaction is present at 
i = 0. However, when dealing with an initially nonempty 
quantum dot we have to assume that the interaction is 
turned on at time t = io avoid initial correlations. 

Green functions in Keldysh formalism — In order to 
describe the time evolution of the system for t > 0, we em- 
ploy the Keldysh formalism.— All single-particle prop- 
erties of interest - such as J and n - can be expressed 
in terms of the retarded and Keldysh component of the 
single-particle dot Green function, 

GZ{t,t') = -ie{t-t')TTpo{d,{t),4{t')}, (9) 
GfAt,t') = -iTTpo[d,it),4{t')]. (10) 

The operators are in the Heisenberg picture with refer- 
ence time to = 0] {...,...} refers to the anticommutator; 
[...,...] denotes the commutator. Finally, the advanced 
Green function 

Gt^^t,t') = ie{t'-t)Trpo{d,{t),dl,{t')} = G]?!{t',t)* 

(11) 

is adjoint to the retarded one. All time arguments are 
positive throughout this paper, t, t' > 0. 

It will prove useful to introduce Green functions that 
are computed w.r.t. three different Hamiltonians: (i) 
the noninteracting, decoupled dot propagator g associ- 
ated with Hq°^ only, (ii) the noninteracting but reser- 
voir dressed dot propagator G^ referring to H with 
iJ™'(t) = 0, and (iii) the interacting and reservoir dressed 
dot propagator G calculated w.r.t. the full H. 

The noninteracting, decoupled dot propagator can be 
obtained by computing dtg{t,t') and dt/g{t,t') and then 
re-integrating with the correct boundary conditions: 

g"'\t,t') = -iQ{t-t')Te-'tt' dt,c{u)^ (j2) 
t') = (t, 0) (1 - 2n)g^''^0, t') . (13) 



T denotes time ordering, and 



(14) 



is the matrix of 'occupancy' of dot states at time t = 0. 

The reservoir dressed but still noninteracting dot prop- 
agator incorporates the presence of reservoirs via appro- 
priate self-energy contributions: 



(15) 



[Ef /i^]^,^(i',t) = T.lLAt')97'^''{t',t)j,At), (16) 



where the noninteracting reservoir propagator gj-^ is 
given in analogy to Eqs. (|12p and p^ . Since we are 
not interested in details of the reservoir band structure, 
we implement a continuous band of infinite width with a 
constant density of states (wide band limit). 



Da{e) = DaC 



-S\e\ 



(17) 



with (5^0^ assuring convergence of the energy integrals. 
This approach is widely used in the literature. Further- 
more, we take the couplings between dot and reservoirs 
to be independent of ka, that is jk^i = lai- This gives 



j:Tit\t) = -iSit'-t)T^it), 
i]^(t',t) = -r„e-^'^°(*'-*)r„^- 



(18) 



1 



SiBh[7TTa{t' -t±i6)]' 

(19) 

where 

Ta.,,,{t) = 7rD,,-/*^AthUt). (20) 
For the derivation of Eq. (jl9|) we exploited the expansion 

1 



1 - 2nfc„ 



■ Ma 



(21) 



where the fermionic Matsubara frequencies uJm are the 
odd multiples of ttTq,, and the series is to be evaluated as 
a principal value for |w„i | oo. Note that the wide band 
limit assumption of a continuous density of states implies 
the limit of infinite reservoir size. Thus, this limit is to be 
performed before further evaluating the dot propagator 
in order to properly define an open system configuration 
in which no recurrence phenomena occur. 
Dyson's equation for reads 



G° = g + G°S]rcs5. 



(22) 



The Green functions and the self-energy in this equation 
are matrices (retarded, advanced, and Keldysh compo- 
nents are ordered in the convention of Ref. l39f). 



(,^0,rot r<O.K \ /vrct vK \ 

^ 1 V - I ^ros ^ros 1 

GO'^'i^' \ ' 



(23) 
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with each block being itself a matrix w.r.t. the dot's 
single-particle quantum numbers. For the retarded com- 
ponent we find 

Go, ret _ ret y^O.rctv^rct ret /'0/l^ 
-9 + ^rcsff 7 l^^J 

where the multiplication abbreviates a summation over 
the dot quantum numbers as well as integration over in- 
ternal times, 

{AB)u\t,t') = Y, dsA,,{t,s)Bj,,{s,t'). (25) 
As EJ^Qg ~ 5{t' - t), the solution is simply 

with Fios = Zq Tq. For the Keldysh component, Dyson's 
equation takes the form 

GO,K _ K /-^O.ret v^ret K /^O.retv'K adv 
- .9 + ^res.9 + "-^ ^res5 

+ G ■ S„s 9 , (27) 

which is solved by 

G^'^{t,t') = -iG°'''"'(t,0)(l-2n)G°''^'^^(0,t') 

^(^0,rct5.K^gO,adv^(^^^,^^ (28) 

The full Green function G finally includes the interac- 
tion. Its components satisfy 

G'°\t,t') = G"^''=*(t,i') + [G"<=ti]r<=tG'0.rctj f^^ ^f^ (29) 
G^{t, t') = -iG'^\t, 0) (1 - 2n)G"'^"(0, t') 

+ [G-'(s^;, + s^)G^^^(^,^'), (30) 

where S is the self-energy associated with iJ'"* (we will 
use the FRG to compute it approximately). The first 
term in Eq. ([30)) governs the decay of the initial dot oc- 
cupancy while the second one describes how a new occu- 
pancy emerges under the influence of the reservoirs and 
the two-particle interaction. 

Physical ohservables — Single-particle properties of 
the system can be expressed in terms of the full Green 
function G. For instance, the (time-dependent) expecta- 
tion value of the occupancy of a dot state easily follows 
from Eq. 

n,{t) = \-'-Gl{t,t). (31) 

Another example is the current of particles leaving reser- 
voir a, 

J^{t) = -iTrpo[H{t),N^{t)], (32) 

where the operators are in the Heisenberg picture. A 
derivation similar to Ref. but in real time space yields 

J„(t) = -Rc Tdt'Tr [!]'„•=*(*, t')G^(i',i) 
Jo 

-G'^\t,t')^l{t\t)]. (33) 

Here, trace and multiplication abbreviate a summation 
over dot quantum numbers only. 



III. FUNCTIONAL RENORMALIZATION 
GROUP 

The functional renormalization group is a quantum 
many-body method which allows to gain insights into the 
physics of interacting fermion and boson systems that 
exhibit a hierarchy of energy scales and/or competing 
instabilities.— In a first step one supplements the nonin- 
tcracting propagation by a flow parameter A, which for 
infrared divergent problems might reasonably be chosen 
as an infrared cutoff. Consequently, one-particle irre- 
ducible vertex functions (effective interactions) acquire 
a cutoff-dependence, and taking the derivative w.r.t. the 
latter yields an exact inflnite hierarchy of flow equations 
- in practice, this can be achieved using a generating 
functional. After truncating the hierarchy - which can 
be done in a strictly controlled way - one obtains a flnite 
closed set of coupled differential equations for the self- 
energy (single- particle vertex), the effective two-particle 
interaction (two-particle vertex) , and possibly higher or- 
der vertex functions (depending on the truncation order). 
The cutoff-free problem is recovered by integrating from 
A = oo where the vertices are known analytically down 
to A = 0. Note that truncation is the only approxima- 
tion; the full dependence on the single-particle quantum 
numbers and times (or frequencies) can be kept. 

Prior applications of the functional RG mainly fo- 
cussed on the regime of linear response where one can 
conveniently resort to the Matsubara formalism (i.e., use 
imaginary times or frequencies). For quantum dot prob- 
lems, two truncation schemes were employed. In the sim- 
plest approximation only the flow of the self-energy E is 
taken into account; E is then frequency independent.— 
Within this truncation all terms to first order in the two- 
particle interaction are kept, while higher-order terms 
are only partially included. We emphasize that the 
RG procedure enhances the quality of the approxima- 
tion beyond perturbation theory (by including contribu- 
tions from an infinite set of Feynman diagrams in a con- 
trolled way; see below). At the end of the RG flow the 
sum of the (frequency-independent) self-energy compo- 
nents merely correspond to an effective noninteracting 
problem,— which allows to gain intuitive insights into 
the observed correlation effects. This static approxima- 
tion can be improved by incorporating the flow of the 
static part of the effective two-particle interaction; see 
Ref. I21I . In a second truncation scheme - which eventu- 
ally includes all second order terms - the full frequency 
and quantum number dependence of the two-particle ver- 
tex and the self-energy are kept, while the flow of the 
generated three-particle vertex is neglected i^^'^^ In this 
improved approximation a significantly larger number of 
coupled differential flow equations needs to be solved; this 
is only possible numerically and requires substantial com- 
putational effort. The second-order scheme was so far 
applied only to the single-impurity Anderson model.— 

The functional RG was recently extended to Keldysh 
frequency space in order to study the steady state of bias 
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voltage driven quantum dots; both approximations men- 
tioned above were employedi^"— The static lowest-order 
truncation is of limited usefulness if one aims at tack- 
ling dots that exhibit Kondo correlations.— One needs to 
employ a higher-order scheme where the self-energy can 
acquire a frequency dependence.— On the other hand, 
the static truncation allows to obtain a comprehensive 
picture of the nonequilibrium steady-state physics of the 
IRLM for small to intermediate interactions: Logarith- 
mically divergent terms that show up in lowest-order 
perturbation theory are resummed consistently, and the 
RG-renormalized tunnel couplings features generic power 
laws with interaction-dependent exponents.— This in 
turn gives rise to highly nontrivial effects such as current- 



voltage characteristics dominated by power laws with 
interaction-dependent exponents. The latter can be com- 
puted to leading order in the interaction.— The simplest 
truncation thus captures hallmarks of correlated charge 
fluctuations within the IRLM. 

Generating functional — In the following we sketch 
how to derive the FRG flow equations for time-dependent 
problems. The procedure is completely analogous to the 
one in equilibrium or steady-state nonequilibrium, and 
more technical details can be found in Refs. Ull, |42, and 
I43I . As a first step, we express the two-time Keldysh 
Green functions via their Keldysh contour functional in- 
tegral representation;^ 



I n ni2PiP2 



(34) 



where ip denote Grassmann fields, and p,p' = ± are the 
usual Keldysh indices referring to the upper and lower 
branch of the Keldysh contour (they can be handled as 
additional quantum numbers for all practical purposes). 
The hat indicates the Green functions before rotation to 
the Keldysh basis (see Sec. |ll]). The interacting part of 
the action is given by 



1,2,1', 2' 



Wl21'2'V'lV'2V'2''01', 



(35) 



using the multi-index 1 = {t,i,p) and the bare two- 
particle vertex 



Ui21'2' =S{ti - t[)d{ti - t2)5{ti - t'2) 



(36) 



where az denotes the z Pauli matrix. The m-particle 
Green function 

Gi...oil'...im' = (~*)™ 

X lT,dfl [t,]... [trn]d\ [t„A . . . 4^^' [ty]) . 

(37) 

can be obtained from the following generating functional 
[the noninteracting part S*^ of the action is given by the 
first term in the exponential of Eq. p4p ]: 



WmAv}) = f-Di^i^ exp{5°-z5'"'-(Vi,77)-('?,^)} 

(38) 

through the derivative 



Gl...in 

(-0' 



5r]i... Srj^ 6r]„ 



-M{v},{v}) 



(39) 



The corresponding functional that generates the one- 
particle irreducible vertex functions is given by the Leg- 
endre transformation 

r({0}, W) = ->v^(W, W)-(^,^)-(^»+(^, [G"]"V) 

(40) 

of the generating functional of the connected Green func- 
tions: 



W^m{v}) = ^n[Wm,{v})]- 



(41) 



Flow equations — If we supplement the free propaga- 
tor by a (for the time being unspecified) flow parameter 
A, i.e., replace G^ -* G^'^, all vertex and Green func- 
tions acquire a A dependence via Eq. (|40p . Taking the 
derivative with respect to A yields the infinite hierarchy 
of FRG flow equations. Their general structure is to re- 
late the A-derivative of the m-particle vertex to a certain 
set of diagrams involving the m + 1-particle vertex and 
lower ones.— This infinite hierarchy can generically only 
be solved by truncating it to a given order. Our approach 
uses the lowest-order scheme for reasons of simplicity; 
more elaborate approximations can in principle be de- 
vised straightforwardly. If the two-particle vertex is set 
to its bare value 



72^(l,2,l',2') = -«Sl21'2' 



(42) 



the only remaining fiow equation is the one for the self- 
energy. It reads 

9A7f(l,l')= E [G^(9A[G°'V)G^UT2^1'2,1',2') 



77=77=0 



■E^2'272^(l,2,l',2'), 
2,2' 



(43) 
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FIG. 1. Diagrammatic representation of the self-energy flow 
equation. The dot indicates the derivative w.r.t. A; the 
slanted line symbolizes the single-scale propagator. Within 
our first-order truncation scheme, the right-hand side is pro- 
portional to a delta function in time, and it is thus sufficient 
to introduce a single time argument t. 



with the so-called 'single-scale propagator' given by 



2,2' 



(44) 



We introduced the star differential operator d'^ which 
acts only on the free Green function not on S"^, 

in the series expansion = G^'^ + G°'^T,^G°'^ + ... . 
The flow equation is depicted diagrammatically in Fig.[T] 
Even though its structure seems Hartrce-Fock-like, it is 
important to stress that our approximation is in no way 
related to any mean-ficldish approach. The latter are 
known to suffer from severe artifacts in low-dimensional 
systems. 

Choice of a cutoff — It is one strength of the func- 
tional RG that it is not bound to a particular cutoff as 
long as the latter fulfills gO'A=~ = o and G°'^=° = G°. 
A crucial step in the extension to noncquilibrium prob- 
lems is to devise a scheme which conserves causality even 
after truncation.— This is guaranteed by the so-called 
hybridization flow; it was successfully used to study the 
steady state of the interacting resonant level model (see 
Rcfs. I23I and [2^ as well as the discussion above). 

The hybridization flow parameter A can be physically 
interpreted as originating from coupling to additional 
auxiliary wide-band reservoirs (one for each dot state) 
with temperature Tfp and chemical potential ^fp.— Af- 
ter rotating to the retarded, advanced and Keldysh basis, 
the corresponding self-energy entering the noninteracting 
but reservoir dressed dot propagator G*^ is determined by 
Eqs. (UHl) and p^: 



I]fp(i',i) = -rfpe-^'^'-(*'-*)A^- 



(45) 



1 



smh[TTTip{t' - t ± i5)] 

(46) 

with Ai'i = ASi'^i. We will frequently suppress the su- 
perscript A (and only reinsert it when crucial to avoid 
misunderstanding). The flow starts at A = 00 where the 
system instantaneously acquires its stationary state (A 
has the physical meaning of a decay rate). The initial 
conditions for the vertex functions are therefore identical 



to those in the stationary state^M^ 7^ °° vanishes for 
n>3, 7^"°° is given by the bare interaction vertex, and 



E:,t^=-(i',i) = ki-oE^ 



(48) 



In prior applications of the hybridization cutoff to the 
IRLM, the temperature of the auxiliary leads was cho- 
sen equal to the physical one. In this work, we employ 
Tfp = 00 since (i) it implies = which simplifies the 
flow equations, and - more importantly - (ii) it avoids 
any imprint of an artificial energy structure from the 
auxiliary leads and might thus retrospectively be a more 
reasonable choice on general grounds. Note that at in- 
finite temperature the chemical potential /ifp no longer 
enters the flow equations. We have checked that steady- 
state results for the IRLM^i are quantitatively unaltered 
if Tfp = 00 instead of Tfp = T is used. 

The Keldysh rotated single scale propagator [whose 
Keldysh component appears on the right-hand side of 
the flow-equation Eq. (|54)) ] is given by 

S' = 9*G= (l + GS])(aAG")(l + EG). (49) 

Since ^aG" = GO(aASfp)G", we find 



with components 



^rct _ ^rct fp ^J^'ct _ _^^ret^ret 

dA 



(50) 



(51) 



as well as 



5^ = G 



■ot "^fp 



dA 



G''°'G^ + zG^G^^^ + G 



G^ + G 



ret fp /-(ad 



dA 



dA 



K 



Gt ^^fp ^ad 



dA 



Finally, the flow equation (|43p translates to 



dA^"^^^ = 0, 



(52) 



(53) 



dAT^'^{t',t) = dAT^f:f{t\t) 

= - E Sf'^{t,t){-tu,,,,,,,.{t))S{t'-t). (54) 

Importantly, the self-energy component E^ does not 
flow. 

Explicit form of the flow equation — In the lowest or- 
der truncation scheme, the expressions for S can be sim- 
plified further: 

T"'\t\t) = S{t' -t)j:''='{t), (55) 
E^(t',i) = 0, (56) 
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and hence 

G''°*(i,i') = -iQ{t-t')Te"' f'*i[f(*i)-»r,cs(ti)-iA+s(ti)] ^ 

(57) 

which in turn entails the multiphcation formula 

G"'\t,t')G'''\t' ,t") = -te{t-t')e{t' -t")G"'\t,t"). 

(58) 

Applying these properties to Eqs. ([STj) and ([52]) eventu- 
ally yields 

S"'\t,t') = {t' -t)G"''{t,t') (59) 

as well as 

S^{t,t') = -{t + t')G^{t,t') 

dhh{G^^\t,h)[{^L + K)G'"''']{h,t') 

+ [G-'(S,t + S|)](t,ti)G^'^(ii,0} 



^ret fp ^ad V 
l9A 



{t,t'), (60) 



where the last term vanishes for our choice Tfp = oo. 

Wc conclude with two comments: (i) The right-hand 
side of the flow equation ([54)) depends on the Green func- 
tion via Eq. (|60|) : thus, the retarded self-energy is fed 
back into the flow through the Dyson equation(s) ([^ 
and (pO)) . This illustrates that Eq. (|54p actually is a dif- 
ferential equation whose integration requires the solution 
of Eq. ([60|) and Dyson's equation for each A. (ii) The self- 
energy originating from the two-particle interaction ap- 
pears in combination with the single-particle part of the 
dot Hamiltonian in all equations at hand. Thus, the for- 
mer can be intuitively interpreted as a time-dependent 
renormalization of the effective single-particle parame- 
ters. Note that this immediately illustrates that the con- 
tinutity equation, which for the most relevant case of two 
leads a = L,R reads 



Mt) + jR{t) = -Y,^,{t) 



(61) 



holds within our approximation. 

Brief summary — Wc have introduced a real-time 
ERG formalism that is capable to describe the time evo- 
lution of interacting quantum dots coupled to Eermi liq- 
uid reservoirs on all scales from transient to asymptotic. 
The dot parameters as well as the level-lead hybridiza- 
tions can carry an explicit time dependence. This equally 
allows to study the transient dynamics of a bias voltage 
driven setup under a time-independent Hamiltonian as 
well as problems that exhibit time-dependent parameter 
variations (e.g., charge pumping). Our derivation of the 
ERG flow equations was kept general: We did not spec- 
ify the number of correlated dot levels and leads or the 
actual geometry determined by the level-lead couplings. 
As an application we will now investigate the relaxation 
dynamics of the IRLM in the presence of a finite bias 
voltage. Results for explicitly time-dependent Hamilto- 
nians will be presented elsewhere.— 




FIG. 2. (Color online) Sketch of the three-site version of the 
interacting resonant level model. 



IV. INTERACTING RESONANT LEVEL 
MODEL 

The dot Hamiltonian of the frequently-mentioned in- 
teracting resonant level model (see Eig.[5]for a sketch) is 
given by 

+ T{d\d2+dld3+U.c.) (62) 

il'"*(0 = 1/(^1712+^2^3), (63) 

where rii = d\di denote three spinless fermionic levels con- 
nected locally through a hopping matrix element r > 
and a Coulomb interaction U (the latter will mainly be 
taken as repulsive; for an exception see Sec. IVH|) . Only 
the central site 2 can be moved in energy by changing the 
'gate voltage' e. The second term in the single-particle 
part of the Hamiltonian is added for mere convenience so 
that e = corresponds to the point of particle-hole sym- 
metry. This term is incorporated into the self-energy; it 
cancels the initial condition Eq. (|47|) of the self-energy 
flow. 

Dot sites 1 and 3 are coupled to left (a = L) and right 
{a = R) noninteracting wide-band leads. Their chemical 
potentials differ by an applied bias voltage V = /^l~/^r ^ 
which we chose symmetrically as /^l = ~Mr = V/2 for 
reasons of simplicity. We focus exclusively on the zero- 
temperature limit. The leads give rise to a self-energy 
term Sjos = + Sr that reads [see Eqs. and p^ ] 



(64) 



-iS{t' - t)Si^i'{6i^i + (5i,3)r, 



where V denotes the principal value, and 



r = r 



Lll 



r 



R33- 



(65) 



(66) 



We have chosen (again for reasons of simplicity) spatially- 
symmetric hoppings and interactions but emphasize that 
these restrictions can be abandoned at the expense of 
minor additional effort. 

In the literature one frequently encounters a field- 
theoretical realization of the IRLM (see Refs. [H, [2^, [sH, 
l46l . and 113 as well as references therein) . The latter con- 
sists of a single spinless fermionic level which is coupled 
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to two leads. A two-particle interaction acts between the 
fermion occupying the level and the ones located at the 
boundaries of the leads. Our microscopic three site dot 
model is equivalent to the field-theoretical IRLM in the 
so-called scaling limit 



T,|e|,|c/|,y«r 



(67) 



where the first and third site of the dot can effectively be 
incorporated into the reservoirs as T is much larger than 
all other (bare) energy scales (see Sec. lVI Al for more de- 
tails). Within our functional RG approach we cannot 
directly treat a single site model with density-density in- 
teraction to the leads. 

The equilibrium physics of the IRLM is dominated by 
an interaction-dependent renormalization of the hopping 
T. First order perturbation theory leads to a logarithmic 
term in the Matsubara self-energy which for e = 
(half filling) reads^i 



M 



u_ 



In 



m 



(68) 



This logarithm is (automatically) resummed by equilib- 
rium functional RG which analytically yields^ii^ 



(t) 

(^) 



-u/{nr)+o{u-') 



|e| « r « r 
T « lei « r 



(69) 



where we have defined r"^ 



^M,A=0 
"12 



T + S 



M,A=0 
23 

The renormalization of the onsite energy e of site 2 is of 
order C7^; the self-energy matrix elements and 
appear in the Green function in combination with T and 
can thus be neglected in the scaling limit. 

The renormalization of the hopping amplitude mani- 
fests physically in the charge susceptibility 



dn2 
Ik 



ttTk 



(70) 



which in turn is governed by the (renormalized) energy 
scale Tk- The latter is the universal scale of the model 
in equilibrium.— We will thus use it as our characteristic 
energy scale. The first-order FRG approximation to Tx 
can easily be computed numerically from ;\;.— 

The nonequilihrium steady- state physics can again be 
solely attributed to a renormalization of the 1, 2 and 2, 3 
matrix elements of the self-energy (with rj^™ + in 
presence of a finite bias voltage V). For |e|,r « V the 
RG flow is cut off by V and hence 



-U I{'itT) + 0{U-) 



(71) 



This leads to a power-law suppression of the current at 
large VM^^^^ 



Jlir (V 



(^) 



-2UI{-nr)+0{U^) 



(72) 



Even richer physics can be observed (i) at e = ±l//2 where 
the differential conductance exhibits a resonance peak, 
and (ii) for left-right asymmetric tunnel couplings to site 
2,29^ In both cases, the bias voltage no longer simply 
acts as an infrared cutoff. We do not consider these pa- 
rameter regimes in the present paper and thus refrain 
from discussing further details. 

After this brief summary of the IRLM's scaling- 
limit physics, we now employ our newly-developed FRG 
scheme to study its real-time evolution. The flow equa- 
tion ([54|) takes the explicit form 



aAE"*(i',t) = 
\ 



-i-5{t' 
2 ^ 



\ 



sUt.ty sUt,t)} 



(73) 



We again emphasize that computing the right-hand side 
of this set of differential equations requires the solution 
of both Eq. ([50)1 as well as the Dyson equation(s) ([29)1 
and We did not succeed in doing this analytically 

(in contrast to the steady-state limit where even the flow 
equation itself was amenable to an analytic treatment^^) 
and thus resort to a numerical solution. In principle^ 
the resulting set of coupled equations can be coded triv- 
ially. In practice, however, the computational effort of 
a straightforward implementation is vast, and no results 
can be obtained in reasonable time. We thus devise an 
efficient algorithm to speed up numerics but do not ap- 
ply further approximations (as does Ref. [28h . The reader 
which is not interested in those computational details 
may skip the next section - only keep in mind that the 
algorithm we are about to introduce is numerically exact. 



V. NUMERICAL ALGORITHM 

To solve the set of equations ([73]), and ((50)) 

numerically we first discretize the continuous time vari- 
able t. The number of time steps M determines the num- 
ber of coupled differential equations that arc eventually 
to be integrated. We are mainly interested in the uni- 
versal scaling limit of the IRLM where the coupling to 
the reservoirs F is much larger than the hopping ampli- 
tude T of two adjacent dot sites. This naturally leads to 
two time regimes characterized by the scales l/F « 1/t. 
Both need to be resolved in order to access all times from 
transient to asymptotic. We thus need to employ small 
discretization steps (« l/F) for small times. However, 
numerical resources limit M to about 10^ (see below), 
so the step size necessarily increases for larger times such 
that eventually the scale 1/t on which the relaxation pro- 
cess towards the steady state takes place can be reached 
(this will be further specified below). 

For discrete time variables t„ (we set to = 0) the full 
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nt-(m+l) 
ri=l 



(74) 



retarded Green function can be decomposed via Eg. (l58l) : results to be numerically exact (typically M ~ 10^ is 

sufBcient}^ numerics can be carried out in reasonable 
time on standard PCs). Thus, no additional approxima- 
tion " such as the often applied^ but in many applica- 
tions uncontrolled generalized Kadanoff-Baym ansata^ 
- is used to integrate the Dyson equation. 

Solving equation (|76|) is particularly simple at J7 = 0. 
Since the self-energy associated with the two-particle in- 
teraction vanishes, there is no need to discretize time. 
One solely has to determine the last term of Eq. ([76|) 
with tn+i and t„ replaced by t and 0, respectively (the 
recursion becomes trivial). The remaining integral can 
be expressed in terms of exponential integrals [the result 
is given by Eq. ()A.7|) with tn+i t, tn ^ as well as A, 



where the product extends over all (discrete) i„ from 
t„j = t to the first time where tm < t' . Thereafter, the 
time arguments of the retarded Green functions differ at 
most by one discretization step, and if the latter is chosen 
small enough, the self-energy can be approximated as 
piecewise constant (which is obviously numerically exact 
in the limit M ->■ oo). This in turn allows evaluation of 
each retarded Green function in Eq. (f74| via Eq. ([57l) : 



and U set to zero]. From G • {t,t) one can calcu- 



G'-^*(t„+i,t„) = -ie-[''o°'+^--^+^f"][*-i-*"], (75) 

where we introduced S = S(S(t), and /ig°* denotes the 
single particle matrix representation of Eq. (|62|) . The 
index i of S''°* indicates the mean value of the retarded 
self-energy in the small interval [i„,i„+i]. The advanced 
Green function G""^^ is related to G"' through Eq. pT|). 

A similar argument applied to the Keldysh component 
yields the recursion relation 

G^ {tn+l,tn+l ) 

= G'^'^^ {tn+l, tn)G^ {tn, tn)G^'^^ (tn, tn+l) 

n-l *yt' *T' 
-iY, / dsi / (is2G'''''(t„+i,si)S^s(si,S2) 



G^ ^{s2,tm+l)G^ ^(<m+l,t,i+l) - H.C. 



^ f f ^'^2G™*(t„+i,Si)E^s(si,S2)C'^^(s2,t„+l), 

(76) 

with the initial condition G^(0, 0) = - 2n). We have 
exploited that G^ and only enter the right-hand side 
of the flow equation ([75]) with equal time arguments. The 
same holds if we want to compute observables: For the 
occupancy ni{t) this is apparent from Eq. ([3T|) : for the 
current Ja{t) of Eq. ([33)l it follows from time locality of 
SJJ^°* in the wide band limit. The integrals in Eq. ([75]) can 
be further evaluated and ultimately expressed in terms of 
exponential integrals. This is discussed in the Appendix 
[see Eqs. (|A.7j) and (jA.lOj) ]. Finally, S^{tn+i,tn+i) can 
be obtained from an analogous procedure. 

The remaining differential (flow) equation (|73|) can 
be implemented straightforwardly using standard Runge- 
Kutta routines. Computing the right-hand side through 
the above recursive procedure scales as M^, and thus 
calculating ni(t) and Ja{t) for a given parameter set ap- 
proximately scales as . We carefully ensure to choose 
the numerical control parameter M large enough for our 



late the noninteracting occupancy ni{t) through Eq. (j31[) 
and the current Ja{t) via Eq. p3|) . Wc will discuss our 
resuhs for C/ = in Sec. IVIAl 



VI. TIME EVOLUTION IN THE IRLM 

A. The noninteracting case 

As an instructive step towards an understanding of 
the nonequilibrium relaxation dynamics of the IRLM we 
study the case U = 0. Throughout this section (and for 
the whole rest of the paper) wc assume that the three 
dot sites are empty at t = but stress that other initial 
states can be considered without any additional effort (as 
long as they fulfill the criteria mentioned in Sec. HI]). We 
will first quantify the heuristic (yet reasonable) statement 
that in the scaling limit defined by Eq. (|57)) our three- 
site model becomes equivalent to the ficld-thcorctical-like 
version of the IRLM where a single site couples to two 
reservoirs. Parameters can obviously be fixed by choos- 
ing 



r 



(77) 



for a given hybridization strength F of the single-site 
model. The latter can be solved analytically at [/ = 0; 
exact expressions for the time-evolution of n and Ja for 
an initially empty and decoupled dot can be found in Ref. 
I32I (see also Ref.|ll|).— In Fig. [3]the occupancy number is 
compared to n2{t) in our three-site dot (the parameters 
are given by r/F = e/F = V/T = 0.025). The two models 
show different behavior at times smaller than l/F (see the 
inset of Fig. [3]), which is the time needed to fill the ini- 
tially empty sites 1 and 3 of the three-site model to their 
steady state values. In the latter the displacement cur- 
rent fills the 1 and 3 site first, whereas in the single-site 
model it gives rise to a nonvanishing first derivative of fi 
for t O^^nSl The difference of the occupancy at small 
times leads to an offset at larger times (which however 
vanishes in the extreme scaling limit F 00). Otherwise 
both models yield identical results. 
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FIG. 3. (Color online) Comparison of the central-site occu- 
pancy between a three-site and a (field-theoretical-like) one- 
site version of the interacting resonant level model for scaling- 
limit parameters e/F = V jV = 0.025, and r/F = 0.025 (the 
hybridization in the single-site model reads F = r^/F; see the 
main text for details). The inset illustrates that both models 
differ only for times ~ 1/F. 
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FIG. 4. (Color online) FRG data for the time evolution of 
the central-site occupancy of the interacting resonant level 
model depicted in Fig. [2l The parameters read t/F = 0.025 
(local coupling), e/Tx = 10 (level position), and V/Tk = 10 
(difference of chemical potentials). The universal equilibrium 
energy scale Tk is renormalized through the local Coulomb 
interaction U. The arrows indicate the steady-state values 
obtained by the nonequilibrium steady-state functional RG 
of Ref . Hi. 



For times t » l/F the occupancy and the current m 
the single-site model are given by (note that Tk = 4r at 



m «n,tat(l + e-^") + -fe 



-2rt 



sin[(e - V/2)t] ^ sin[( 



{e-V/2rt (e 



+ V/2)H ) 



Ja{t) 



^1 - 2n{t) H — I - arctan 



2f 



(78) 



(79) 



where we have introduced the stationary occupancy 
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FIG. 5. (Color online) The same as in Fig.|4]but now showing 
the current leaving the left reservoir. The negative current 
found at small times was already observed in Refs.[2^ and |3^. 



J^stat = J + ^ E arctan ( j . (80) 

The corresponding expressions of the three-site model at 
large F are identical (if n and F are replaced by n2 and 
T^/F, respectively). Equation ([78|) illustrates that the 
long-time behavior of the occupancy exhibits two oscil- 
latory terms with frequencies e ± V /2 as well as an ex- 
ponential relaxation with rates 4F and 2F; the latter is 
accompanied by a power-law correction 1/t. The cur- 
rent characteristics are similar. For later reference we 
note that in the limit |e ± V/2\ » Tk the current Jq, 
only shows the single frequency e - fJ-a since the second 
frequency (originating from the occupancy term only) is 
suppressed. 



B. The interacting case 

After these prerequisites we now turn to the interact- 
ing case. Our objective is twofold. We first discuss our 
results from a pure physical perspective, assuming that 
our FRG scheme yields accurate data as long as U is 
not too large. This is certainly reasonable in light of the 
method's success to capture correlation effects within the 
IRLM in equilibrium and steady-state nonequilibriumi^ 
Thereafter, we relate our findings to real-time RG calcu- 
lations (Sec. IVI C|) . The RTRG access to nonequilibrium 
dynamics differs strongly from the FRG approach, and so 
do the applied (controlled) approximations. Thus, com- 
paring with RTRG data provides a highly nontrivial test 
for the newly-developed FRG framework. 

Our discussion of the interacting case focusscs exclu- 
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sively on the limit |e ± V I2\ » Tr-, which is a posteriori 
motivated by the fact that the physics can be interpreted 
in a simple way. Moreover, the far-from-equilibrium 
case V » Tk is undoubtedly most intriguing on general 
grounds. We again point out that the FRG is not bound 
to any parameter regime - it can flexibly describe the 
whole crossover from V « Tk to V » Tk (to V »T). 

The results for the time dependence of the occupancy 
n2{t) and the current Jl(0 at different U but fixed t/F = 
0.025, e/Tk = 10 = V/Tk are depicted in Figs. H and 
[SJ The axes are scaled using the universal equilibrium 
energy Tk ■ Note that e and V need to be varied with U in 
order to keep e/Tk and V/Tk fixed. Whereas the steady- 
state values of n2 and Jl are nonuniversal, the oscillation 
frequencies depend only very weakly on the interaction 
strength U (the relative position of maxima and minima 
remain unaltered). In analogy to the noninteracting case, 
only a single frequency |e + V/2\ governs the current in 
the limit |e± V^/2| » Tk while both |e±F/2| are manifest 
in the occupancy. 

It is a distinct advantage of the first-order FRG frame- 
work that it allows to gain intuitive physical insights from 
investigating the (time dependence of the) rcnormalizcd 
effective single particle parameters. Figure |6] shows the 
renormalized onsite energy e'^™(t)-e=I]^'2*'''="eR (note 
the double logarithmic scale). For times larger than l/F, 
it does not flow to leading order in U. The renormal- 
ized hopping amplitudes rj;™ - t = S'^j''^"^ e C and 
T2'§" - T = ±lf^^° € C are depicted in Fig. [7] (recaU that 
rj;™ = only in absence of a bias voltage) . They begin 
to oscillate very quickly around their steady-state val- 
ues. The renormalized hopping between sites 1 and 2 
(between 2 and 3) only exhibits the smaller frequency 
e-V/2 (the larger e + V/2), which was checked via Fourier 
transformation. Thus, the oscillation frequencies e±V/2 
which govern observables at J7 = are not altered through 
the Coulomb interaction; the latter merely leads to U- 
dependent phase shifts. 

This observation indicates that the time evolution of 
n2 and Ja in the interacting case might be described by 
an analytic expression similar to the one at J7 = if the 
latter is supplemented by ?7-dependent prefactors, rates, 
and phase shifts. For a further investigation we make the 
ansatz (sticking to t» l/F for simplicity) 

n2{t) = ao + + a^e-'"''^-^^ " ^ ^ 

(e-"-|:) t 

+ 056 2 ■ 

(81) 

This is precisely the form of Eq. (|78p if it is generalized to 
account for left-right asymmetric bare hoppings (ri2 i= T23 
is required at > 0); one can easily show that at J7 = 0, 
the coefficients are then given by uq = ai = n2.stat{U = 0), 
a2 = 04 = + r23, 03 = 2ri2/7r, and 05 = 2r23/7r, where 
Tij = T^JT. At [/ > 0, the Oi serve as fit parameters. Our 



T^t=2T^/r 




FIG. 6. (Color online) Time dependence of the renormalized 
onsite energy e'^™(t)-e e R (of the central site 2). The parame- 
ters are the same as in Fig.H t/F = 0.025, e/Tx = V/Tk = 10, 
and U/T = 0.1. Note the double logarithmic scale. Kinks in 
the graph are due to the time discretization. 
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FIG. 7. (Color online) The same as in Fig. [6] but for the 
renormalized hoppings r^"' - r = E'^j' '^" C between sites 1 
and 2, and r23° - r = Ejs* '^" ^ C between sites 2 and 3. 



conjecture that a crucial part of the interaction effects 
can indeed be incorporated by taking the noninteract- 
ing functional form but renormalized parameters is then 
strongly supported by the fact that (i) the fitting error 
does not increase with U (sec Fig. [8] for an illustration), 
and (ii) the interpretation of the coefficients is consistent 
with the noninteracting case: oq and ai are very close 
to the [/-dependent steady-state occupancy as directly 
obtained by steady-state functional RG^ or averaging 
the real-time data for large times; likewise, the rates 02,4 
and the asymmetry ratio a^/a^ agree with the asymp- 
totes of the renormalized F12, F2,3, and Fi2/r2,3 (see Tab. 
|l|. In passing, we note that the phase cf> increases from 
(/)(C//F = 0) = to ^(J7/F = 0.2) ^ 0.26 for the parameters 
of Fig. m and thus shows a significant [/-dependence. 
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FIG. 8. Relative difference of the data for the occupancy 
shown in Fig. 2] and a fit to Eq. (|81|l . Small deviations even 
in the noninteracting case originate from F being finite and t 
being only roughly an order of magnitude larger than Tk- 



FIG. 9. (Color online) Comparison of FRG and real-time 
RG data for the occupancy n2(t) at r/F = 0.025, and e/Tx = 
V/Tk = 10 (the parameters of Fig.©. The curves at U/T = 0.1 
are almost indistinguishable. 



0.02 



C. Comparison to real-time RG 

In this section we relate the data from our newly- 
developed real-time FRG scheme to independently ob- 
tained results. As a first consistency check we com- 
pare the long-time asymptotes of the current and the 
occupancy with a direct calculation using lowest-order 
nonequilibrium steady-state functional RG.— As illus- 
trated in Figs, m and [5] - where arrows indicate values 
extracted from the independent steady-state framework 
- the agreement is quantitative. 

The time evolution of the IRLM exposed to a bias volt- 
age was recently studied using real-time RG.— The 
RTRG is based on approximations that are controlled at 
small U and t. In equilibrium and steady-state nonequi- 
librium FRG and RTRG data was shown to agree on a 
quantitative level at small [/.— For the discussion below 
one should bear in mind that it is not straightforward 
within the RTRG to obtain prefactors correctly to or- 
der U due to technical (truncation) subtleties.— Rates 
and exponents of possible power-law corrections, how- 
ever, come out correctly to order U This is different 
in the truncated FRG used here where it is guaranteed 
that the prefactors are correct to leading order in the 
interaction. 



U/r (Fi2 + F23)/rA' Fi2/r23 a2/TK aijTK as/as 

0.500 1.000 0.504g:|05 0.496o:l93 1-056 

0.1 0.453 1.073 0.463g:li8 O.US^iH 1.071 

0.3898:if 0.405J!:gi 1.167 



0.2 



0.417 



1.149 



TABLE I. The fitting parameters 02/4 in comparison to the 



renormalized relaxation rates Tij - — p 
error is with respect to a 68% confidence level 
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FIG. 10. (Color online) The same as in Fig. |9] but for the 
current. 



A comparison of the occupancy and the current ob- 
tained by FRG/RTRG is shown in Figs. |9] and [lOl the 
agreement is satisfying. Small differences may be at- 
tributed to the fact that our parameter set is not in 
the extreme scaling limit (for a quantitative analysis of 
this, see Sec. lVI A)) , while the real-time RG is directly set 
up in this regime. Deviations of similar magnitude were 
observed in prior steady-state calculations)^ Moreover, 
RTRG prefactors are not controlled to leading order, and 
this issue becomes more relevant at larger interactions. 

It is a key advantage of the real-time RG that one 
can derive approximate analytical expressions for the 
long time behavior of the occupancy and the current in 
the regime |e ± V/2\ » Tk even in presence of (not too 
large) interactions.— 1^ In the unbiased case, the current 
at sufficiently large t reads (observe that V = implies 
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J{t ^ oo) ^ 0) 



T, 



ZTT 



-r,t/2 



cos(ei) 



It 



(82) 



where T, T^, and e denote the RTRG-renormahzed de- 
cay rates and level position. One particularly finds an 
interaction-dependent power-law correction with 
g = 2U/Tr + 0{U^) [compare to Eq. (HH) for [/ = 0]. 
Real-time RG predicts power-law corrections also for 

V > om^ 

It is an important test for our approximate functional 
RG approach - which is based on a truncation scheme 
guided by perturbation theory in orders of f7 - to inves- 
tigate whether the RG resummation not only leads to an 
exponential decay in time with [/-dependent decay rates 
(as was already proven above) but whether it simulta- 
neously gives rise to power laws in t with [/-dependent 
exponents. The predicted power-law corrections to some 
of the exponential decay term o^^i'^^ are subleading correc- 
tions. Achieving a rather good fit of our functional RG 
data with the ansatz of Eq. (|8T|) thus does not exclude 
the presence of such terms. 

We now analyze our FRG data for potential power-law 
contributions. For reasons of simplicity, we focus on the 
unbiased case where the hopping amplitudes renormal- 
ize symmetrically. Note that due to the structure of Eq. 
(|82|) a fitting procedure can easily shift the influence of 
a power law to a small variation of the relaxation rates. 
We thus do not fit the relaxation rates but rather take 
them as their steady-state renormalizcd values (see Sec. 
IVIBl in particular Fig. [7] and Tab. Likewise, we re- 
place e by our steady-state renormalizcd level position 
(see Fig. in]). Thereafter, we fit the FRG data to the form 
predicted by real-time RG [Eq. ([82)) ] employing oi and 
g as fitting parameters. We emphasize that the prefac- 
tor Tfc of the oscillatory term in Eq. (|82|) does not suffer 
from the above mentioned subtleties of the RTRG, but 
«! does and must thus be fitted. As illustrated in Fig. 
m] the [/-dependence of the extracted exponent g agrees 
with the real-time RG prediction within the error bars. 
This indicates that our functional RG indeed achieves 
a resummation which simultaneously leads to an expo- 
nential time decay as well as power-law corrections, both 
featuring U -dependent exponents. Verifying this analyti- 
cally (e.g., using the cutoff introduced in Ref. [13) will be 
subject of a future investigation. 
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FIG. 11. Interaction dependence of the power-law exponent 
extracted from our FRG data by fitting to the RTRG predic- 
tion of Eq. (|82p at e/T/f = 30 (see the main text for details). 
The solid line shows the leading order term 2U /vr. The error 
bars refer to a 68% confidence level. 



improved noninteracting blip approximation.— Relating 
our data to those results provides an additional test for 
the newly-developed time-dependent functional RG ap- 
proach. 

The spin-boson model describes two fermionic states 
tunnel-coupled by A and separated by an energy E which 
interact with a bath of bosons: 



E 

HsB = —<^z 



where aip are bosonic lowering (raising) operators, a^jz 
denote the Pauli matrices, gq describes the momentum 
q dependent coupling to the phonon bath, and the 
dispersion of the latter. For the so-called ohmic case 
where the spectral density reads 



A 



^T.9q{aq+al), (83) 



J(w) = TT ^ gp{uj -ujq) = 27raa;e""^'^ 



(84) 



with ujc being a high-frequency cutoff, one can map the 
spin-boson model to the unbiased {V = 0) IRLM.^- Tj^g 
parameters are related through 



e = E 
U/inT) = (l-2a)/2 



(85) 
(86) 



VII. RESULTS FOR THE OHMIC SPIN-BOSON 
MODEL 

As a second application we focus on a special param- 
eter regime for which the IRLM can be mapped onto 
the so-called ohmic spin-boson modelj^ The latter con- 
stitutes one of the basic models used to study decoher- 
ence and relaxation phenomena in quantum systems, and 
its relaxation dynamics was previously investigated by a 
variety of methods, in particular field theor}*^ and an 



8uJc 



(87) 



The occupancy of the IRLM determines the expectation 
value of cr, via 



{'yz) + i 



where we can identify n{t) = n2{t) of our three site dot 
model in the scaling limit. We again assume an initially 
empty dot, i.e. a spin pointing in -z-direction at t = 0. 
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For e = E = the ohmic spin-boson model suppos- 
edly exhibits a coherent-incoherent transition at a = 1/2 
(where the interaction U of the IRLM changes from re- 
pulsive to attractive).— Field theory^ and an improved 
noninteracting blib approximation^ both predict the re- 
laxation rate and the frequency Q (in the coherent phase) 
to be related through 



4r 



(89) 



We have dropped second-order terms which are not con- 
sistently included within our truncation scheme. 

FRG results for the time evolution of |(cr2(i))| = 
\2n2{t) - 1| are shown in Fig. [121 As expected, the ex- 
ponentially damped oscillations (indicated by the dips) 
present for U > disappear at f7 < as the dynamics turn 
incoherent. We scale the y-axis logarithmically (i) to em- 
phasize the exponential damping, and (ii) to illustrate 
the coherent-incoherent transition at U = more clearly. 
Note that only a single oscillation period manifests even 
for the largest interaction. This is qualitatively consis- 
tent with Eq. (|89|) which postulates a frequency AUT/T; 
in the scaling limit, this is a very small number if U is 
not too large. For a quantitative analysis, we fit our FRG 
data to the predicted long-time behavior: 



fi2(t) = 0.5 + ae"^*cos(f7i) 



(90) 



The frequency extracted from the FRG framework 
agrees with Eq. ((M)) even for fairly large values of U (see 
Fig. [T5)) . One should note that the coherent-decoherent 
transition is exclusively driven by the interaction (i.e., by 
the coupling of the spin to the bosons). Reproducing the 
transition as well as the relation between the relaxation 
rate and oscillation frequency (the latter being propor- 
tional to Ul) in the coherent regime constitutes another 
stringent test of our approximate approach. 



VIII. CONCLUSION 

We developed a functional RG approach to study time- 
dependent electron transport through quantum dots with 
a few correlated degrees of freedom coupled to Fermi liq- 
uid reservoirs. It allows to investigate the relaxation dy- 
namics out of an initial nonequilibrium state into the 
steady state driven by a bias voltage. Additionally, the 
dynamics of Hamiltonians with time-dependent param- 
eters can be tackled. The hybridization flow parame- 
ter was used, and the hierarchy of functional RG flow- 
equations was truncated to lowest order. 

We applied this approach to the interacting resonant 
level model with time independent parameters and com- 
puted the time dependence of the occupancy of the dot 
level and the current through it. We devised an effi- 
cient algorithm to solve the Dyson equation as a key step 
for a numerically exact implementation of the flow equa- 
tions. The relaxation dynamics of the system exposed 
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FIG. 12. (Color online) FRG calculation of the expectation 
value I (a^it)) \ = \2n2{t) - 1| for the ohmic spin-boson model. 
In absence of a bias voltage, the latter can be mapped to the 
IRLM, and our fermionic FRG scheme is directly applicable. 
In terms of the IRLM, the system parameters read r/F = 0.025 
and e = V = 0. Note that the y-axis is scaled logarithmically 
and that the times are given with respect to F instead of Tk- 




m 0.4 

u/r 

FIG. 13. (Color online) Oscillation frequency governing the 
relaxation dynamics of the spin-boson model as a function of 
U for r/F = 0.025. The functional RG result is compared to 
the prediction of Eq. (|89|) . The error bars result from a 68% 
confidence level fitting of the frequency. 



to a bias voltage is dominated by an exponential decay 
with two different rates and oscillations with two differ- 
ent frequencies. While the decay rates are significantly 
renormalized by the interaction, the frequencies are al- 
most unaffected. In addition, we observed a power-law 
correction with an interaction-dependent exponent. Re- 
markably, our functional RG procedure thus (automati- 
cally and consistently) carries out two different types of 
resummations: one leads to the combined appearance of 
the time t and interaction U in the argument of expo- 
nential functions; the other gives rise to a power-law in 
t with a t/-depcndent exponent. This is highly nontriv- 
ial. We compared our results to those obtained within a 
recently-developed approximate real-time RG approach 



15 



Resi 



11 



22 



12 



13 



1 1 + 



(t>Jl-lJ2)("l-'^3) 



)([iJ3-lJ2) "3-"2 "3-'^2 (l^3-"l)("3-^2) 



33 



21 



31 ReSl3,„(ri'^2,T23 ^ (T12)*, 



23 



32 



TABLE II. Residues of Eq. ^K^. 



ReSl2,n(Ti'^ -* T^i) 



Res23,n(r2l, - (ral,)*) 



and observed good agreement. 

We then exploited a mapping of the unbiased IRLM 
to the ohmic spin-boson model. The latter features two 
localized states; if their energy is equal, one expects a 
coherent-dccohcrcnt transition, which is confirmed by our 
calculation. The FRG prediction for the relation between 
the relaxation rate and the oscillation frequency agrees 
with those obtained by a field theoretical approach and 
an improved nonintcracting blib approximation. Note 
that this relation is a manifestation of strong-coupling 
physics: The explicit scale set by the energy difference of 
the two levels is zero, and the coherent-dccohcrcnt tran- 
sition is purely driven by the interaction (the oscillation 
frequency is proportional to U). Reproducing these re- 
sults hence provides another stringent test for our newly- 
developed method. 

The time-dependent functional RG approach directly 
allows to tackle explicitly time-dependent Hamiltonians. 
This constitutes a highly-active field of current research, 
and one can readily envisage a vast number of appli- 
cations. From a theoretical perspective, it is intriguing 
to investigate the time evolution of systems with corre- 
lated initial density matrices towards their steady-state 
through the following protocol: one starts with an un- 
correlated initial density matrix; the system then relaxes 
towards a steady-state which contains correlations; there- 
after, the system is quenched, and one studies the relax- 
ation process out of the now correlated initial state. From 
a practical point of view, one can investigate correlation 
effects on quantum pumps (which require a periodic vari- 
ation of the dot parameters). 
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APPENDIX 



To solve the recursion Eq. (1761) wc need to compute 
the integral 



j dsi J (iS2G'°'(tn+i,Si)S^s(si,S2)G'"^^(s2,im+l) 

(A.l) 

with I]^g(si,S2) as in Eq. (|55|) . The exponential of the 
matrix in Eq. (|75p can be evaluated: 



Glf{t,t') = -z^ Res,,, e-'"' (A.2) 
1=1 

with ReSij,; and loi being the residues and the poles of 



(A.3) 



Introducing the effective parameters 



A _ ^ , v^rct 
A 



T,, = T + i:r 



(A.4) 



~23 - T + ^t,23 

allows to express the poles as 



t,12i 
ret 



wi = e'^-i(r + A), c^2/3=^(f 



'A 



iV - 2iA T 



i,A\2 



)2+4|^A|2+4|^A| 



(A.5) 



The corresponding residues are given in Tab. |TT1 



16 



To compute the integral Eq. (jA.ip one substitutes T = ti +t2 and At = t2- ti 



J dti J dt2 — > - J dAt J dT+ J dAt J dT + J dAt J dT 

tm-tn + l 2t„,-At tm + l-t„ + l 2t,„-At tm-tn 2t„+At 



(A.6) 



for i„+i - tn > tm+i - tm (thc Opposite case follows analogously). Because of the principal value involved it proves 
advantageous to separate the problem into the two cases n = m and n + m. The case n = m also includes the solution 
of the noninteracting system, where one can choose a single discretization step; we will illustrate it first. With the 
above substitution one can write the i,j matrix element of Eq. (|A.1[) in terms of exponential integrals)^ 



tn+l 



J dsi J ds2G"*(t„+i,Si)EK^(si,S2)G'^''"(s2,t„+l) 



t„ t. 

-5 



= hm - ^ Res,„,„Res;„ e-^---*-! 

ij r.,™.l,2,3 



2t„+i+At t„+i-t„ 2t„+i-At 

J dAt J '^^^ f '^^^ f '^^ 

tn-tn+i 2t„-At l5 2t„+At 

op 

V Res- Res* p-»(""*"+i-"„*"+i) 



At 



1=1,2,3 



iAWnrnTT 

where one has exploited that 



(A.7) 



- e 



- e 



lim Ei{-x5) -Ei{yd) = -log(-a;) +\og{y) , 



and defined 



(A., 



(A.9) 



In the indices of the residues in Eq. (|A.7p one has to replace a = L by 1 and a = i? by 3. For the case m + n one 
analogously finds 



J dsi J dS2G"*(^„+l,Sl)S^;,(Sl,S2)G^^"(s2,Wl) 



1 



V] R,eSjQ TjHes,,-^ , 



tm+i-t„+i 2t„+i+At fm-t„ 2t„+i-At tm+i-t„ 2t„,+i-At 

j dAt j dT+ J dAt J dT+ J dAt J dT 

tm-tn + 1 2t„-Af t,„ + i-f„ + i 2t™-Af 

■( t * t ) 2-^ 

R,esiQ,_^Hes^-Q 



2f„+At 



AfgiiTAw„,„g-iiAt(w„+w*J 

At 



1-1,2,3 



+ e*'^"""*"i;i(-i[/,t„ - w„][t™ - 1„]) - e''^'""'"*'"i;i(-i[/i„ - w„][t™ - t„+i]) 



(A.IO) 
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